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ABSTRACT 



We have performed a 70 billion dark-matter particles N-body simulation in a 2 /z~^ Gpc periodic box, using the concordance, cosmo- 
logical model as favored by the latest WMAP3 results. We have computed a full- sky convergence map with a resolution of A^ ^ 0.74 
arcmin^, spanning 4 orders of magnitude in angular dynamical range. Using various high-order statistics on a realistic cut sky, we 
have characterized the transition from the linear to the nonlinear regime at ^ ^ 1000 and shown that realistic galactic masking affects 
high-order moments only below £ < 200. Each domain (Gaussian and non-Gaussian) spans 2 decades in angular scale. This map 
is therefore an ideal tool for testing map-making algorithms on the sphere. As a first step in addressing the full map reconstruction 
problem, we have benchmarked in this paper two denoising methods: 1) Wiener filtering applied to the Spherical Harmonics decom- 
position of the map and 2) a new method, called MRLens, based on the modification of the Maximum Entropy Method on a Wavelet 
decomposition. While the latter is optimal on large spatial scales, where the signal is Gaussian, MRLens outperforms the Wiener 
method on small spatial scales, where the signal is highly non-Gaussian. The simulated full-sky convergence map is freely available 
to the community to help the development of new map-making algorithms dedicated to the next generation of weak-lensing surveys. 

Key words. Methods: N-body simulations; Cosmology: observations; Techniques: image processing 
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1. Introduction 

Weak gravitational lensing, or "cosmic shear", provides a 
unique tool for mapping t he matter density distribution in the 
Uni verse (for rev i ews, s ee iRefregieil (l2003l) . iHoekstral (l2003l) 
I and iMunshi et alJ (l2006h ). Current weak-lensing surveys cover 
. altogether about 100 square degrees and have been used to 
■ measure the amplitude of the matter power sp e ctruni and other 
] cosmological parameters (see iBenjamin et alJ (l2QQ7l) . iFu et alJ 
' (1200 8) and references therein). A number of new instruments 
' are being planned to carry out these surveys over wider sky 
! areas (PanSTARRS, DES, SNAP and LSSTfl or even over 
' the full extragalactic sky (DUN^. These wide-field surveys 
will yield cosmic-shear measurements on both large scales, 
where gravitational dynamics is in the linear regime, and small 
scales, where the dynamics is highly nonlinear. The compar- 
ison of these measurements with theoretical predictions of 
the density field evolution will place strong constraints on 
cosn iological pa rameters, including dark energy parameters 
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(e.g.lHu & Tegmark( 1999h.lHutereill2002b . lAmara & Refregied 
(l20Q6h and lAlbrecht & BernsteinI (l2007h ). On small scales, the 
highly nonlinear nature of the density field ensures that pre- 
dictions based on analytic calculations are prohibitively difi^i- 
cult and requires the use of numerical simulations. N-body sim- 
ulations have thus been used to simulate weak-lensing maps 
across small patches of the sky, using the flat sky approximation 
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e.g. iJain et al.l (l2000h . iHamana et al.l (l2Q0Th and lWhite & Vald 
|2004|) ). The simulation of full-sky maps in preparation for fu- 
ture surveys involve a wide range of both mass and length scales 
and is challenging for current N-body simulations. The range of 
scales involved also requires the development of efiftcient algo- 
rithms for deriving a mass map from true noisy data sets. These 
algorithms need to be well-suited to both the large-scale sig- 
nal, which is essentially a Gaussian random field, and those on 
small-scales, where it is highly non-Gaussian and exhibits lo- 
calized features. 

In this paper, we used a high resolution N-body simulation 
to construct a full- sky weak-lensing map and test a new map- 
reconstruction method based on a multi-resolution technique. 
For this purpose, we use the Horizon simulation, a 70 billion 
particle N-body simulation, featuring more than 140 billion cell 
in the AMR grid of the RAMSES code (Teyssier 2002). The 
simulation covers a sufi&ciently large volume (L^ox = 2h~^Gpc) 
to compute a full- sky convergence map, while resolving Milky- 
Way size halos with more than 100 particles, and exploring small 
scales into the nonlinear regime (see Sect. [2]). This unprece- 
dented computational eflTort allows us, for the first time, to close 
the gap between scales close to the cosmological horizon and 
scales deep inside virialized dark-m atter haloes. A similar eflfort 
at lower resolution was presented by lFosalbaet al.l (l2QQ8l) . 

The dark-matter distribution in the simulation was integrated 
in a light cone to a redshift of 1 , around an observer located at 
the centre of the simulation box (see Sect. [3]). This light cone was 
then used to calculate the corresponding full- sky lensing conver- 
gence field, which we mapped using the Healpi}0 pixelisation 
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Fig. 1. Full- sky simulated convergence map derived from the Horizon Simulation. Its resolution of 200 million pixels has been 
downgraded to fit the page. The various inserts display a zoom sequence into smaller and smaller areas of the sky. The pixel size is 
0.74 arcmin^. 



scheme (iGorski et al.ll2005h with a pixel resolution of A^ ^ 0.74 
arcmin^ (^side = 4096), and added "instrumental" noise for a typ- 
ical all-sky survey with 40 galaxi es per arcmin^, as expe cted for 
example for the DUNE mission (JRefregier et al.l 120061) . Using 
an Undecimated Isotropic Wavelet Decomposition of this realis- 
tic simulated weak-lensing map on the sphere, we analyzed the 
statistics of each wavelet plane using second, third and fourth 
order moments estimator (Sect.©. We then applied, in Sect.[5j 
a multi-resolution algorithm to filter a fictitious simulated k data 
set based on an extension of the wavelet filtering technique of 
IStarck et aP (2006b). We characterised the quality of the recon- 
struction using the power spectrum of the error map and com- 
pare this to the result of standard Wiener filtering on the sphere. 
Our results, summarised in Sect. [6l illustrate the virtue of high 
resolution simulations such as the one reported here to prepare 
for future weak-lensing surveys and to design new map-making 
techniques. 



2. The Horizon N-Body simulation 

This large N-bo dy simul a tion was carried out using the 
RAMSES code (iTevssied l2002h for two months on the 
6144 Itanium2 processors of the CEA supercomputer BULL 
Novascale 3045 hosted in France by CCRT0. RAMSES is a 
parallel hydro and N-body code based on the Adaptive Mesh 
Refinement (AM R) techniques. Using a parallel version of the 
grafic package (JB ertschin ge j |200 lb . we generated the initial 
displacement field on a 4096^ grid for t he cosmological pa - 
rameters from the WMAP 3rd year results (ISpergel et al.ll2007l) . 
namely Q^ = 0.24, Qa = 0.76, Q.^ = 0. 042, n = 0.958, Hp = 73 
km/s/Mpc and erg = 0.77. We used the lEisenstein & Hul (Il999h 
transfer function, which includes baryon oscillations. The box 
size was set to 2 Gpc/h, which corresponds roughly to the co- 
moving distance to an object at z - 0.8. We used 68.7 billion 
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particles to simuate the dark-matter density field, yielding a par- 
ticle mass of 7.7 x 10^ M© and 130 particles per Milky Way halo. 
This large particle distribution was split across 6144 individual 
files, one for each processor, a ccording to the RAM SES code do- 
main decomposition strategy (iPrunet et al.ll2008l) . Starting with 
a base (or coarse) grid of 4096^ grid points, AMR cells are re- 
cursively refined if the number of particles in the cell exceeds 
40. In this way, the number of particles per cell varied between 
5 and 40, so that the particle shot noise remained at an accept- 
able level. At the end of the simulation, we had reached 6 levels 
of refinement with a total of 140 billion AMR cells. This cor- 
responds to a formal resolution of 262 144^ or 7.6 /z"^ kpc co- 
moving spatial resolution. Parallel computing is perfomed using 
the MPI message-passing library, with a domain decomposition 
based on the Peano-Hilbert space-filling curve. The work and 
memory load was adjusted dynamically by reshuffling particles 
and grid points from each processor to its neighbors. The sim- 
ulation required 737 main (or coarse) time steps and more than 
10^ fine time steps for completion. 

3. Light cone and convergence map 

Born's unperturbed-trajectory assumption for all neighboring 
light rays is a good approximation in the linear regime of 
structure formation, but is inaccurate in the nonlinear regime. 
Consequently, distortion eff'ects of lensing beyond the first order 
canno t be simulated reliably. As shown by IVan Waerbeke et all 
(I2OOII) . the Born approximation also introduces a relative er- 
ror in the skewness of the signal of aproximately 10% on large 
scales where the convergence is Gaussian, and about 1% on 
small scales in the nonlinear regime. We therefore implemented 
a multiple-lens ray-tracing method that can be applied more gen- 
erally than Born's approximation. 

We constructed a light cone by recording, at each main time 
step, the positions of particles within the boundaries of a pho- 
ton plane: this plane moved at the speed of light towards an ob- 
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Sky mask 




Fig. 2. Map of the cut-sky used in Sect.|4]to compute high-order 
moments. 



server, who was located at the centre of the box. Our method 
was developed from the one presented by lHamanaet al.l (l2QQli) . 
This method produced 348 slices in the light cone, spanning the 
redshift range [0,1]. Due to the large size of the simulated vol- 
ume, the effect of periodic replications of the computational box 
are minimized. Each slice was then transformed into a full-sky 
Healpix map (nside = 4096) of the average overdensity using 
a simple "Nearest Grid Point" (NGP) mass projection scheme. 
The density slices thus represented our physical model of the 
lens screens used in the ray-tracing procedure. We note that there 
is no unique procedure for generating a band-limited harmonic 
representation of each slice of particles. We chooe to use an NGP 
interpolation because it is a good compromise between filtering 
and aliasing, and remains localised in configuration space. More 
sophisticated interpolation schemes have been developed in the 
context of either 3D par ticle distribiitions (IColombi et al.ll2008|) 
or 2D continuous fields (JBasak et al.ll2008l) . which, however, re- 
main impractical in significantly large simulations. 

After an interpolation kernel has been chosen, all fields (lens- 
ing potentials and displacement fields) are computed from the 
NGP interpolation mass slices at each redshift using a spher- 
ical harmonic decomposition. The resampling of the displace- 
ment fields outside the pixel centres (as required in a multi-lens 
method) is completed using a local linear-interpolation scheme 
(using covariant, second derivatives of the potential); this last 
interpolation has the same spectral behavior (and thus the same 
aliasing contamination) as the NGP-interpolated mass slices, 
and we therefore do not need to use a higher-order resampling 
scheme (since the calculation of the potential requires two sets 
of integration over the mass distribution, while the interpolation 
of the displacement field corresponds to a second-or der deriva- 
tive). We provi de mo r e deta ils in Appendix |Al (see 'Jai n et al.l 
(I2000h . Ha mana et al.1 (l200ll) and White & Vale (2004) for al- 
ternative approaches). We assumed that the background galaxies 
are within a single source plane located at redshift z^ = 1 . The fi- 
nal convergence map was computed using our multiple-lens ray- 
tracing scheme, for which spherical geometry pre cludes the use 
of small angle approximations (as in lDas & Bode (2008) ) espe- 
cially in the neighborhood of the poles; full rotation matrices for 
each light ray must therefore be computed from the displacement 
fields at each redshift. 

The resulting full- sky Healpix map with a pixel size of A^ ^ 
0.74 arcmin is shown in Fig.[T] with small inserts to highlight the 
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Fig. 3. Moments of the convergence as a function of the average 
multipole moment on each wavelet scale. The variance, skew- 
ness, and kurtosis are shown as black, blue, and red lines, re- 
spectively. Solid lines with error bars corresponds ro a full-sky 
analysis, while dotted lines correspond to our cut-sky analysis. 



large dynamical range achievecQ. The particle shot noise corre- 
sponding to our 70 billion particle run has a small impact on the 
map. As shown in Fig. IH the particle shot noise is well below 
the expected instrumental noise, and even sufliiciently low to be 
ignored in the spectral analysis of the signal. 



4. High-order moments and realistic sl<y cut 

In Fig. [H the signal appears as a typical Gaussian random field 
on large scales, similar to the Cosmic Microwave Background 
map seen by the WMAP satellite ( Spergel et al. 2007). On small 
scales, the signal is clearly dominated by clumpy structures (dark 
matter halos) and is therefore highly non-Gaussian. To character- 
ize this quantitatively, we performed a wavelet decomposition of 
our map using the Undecimated Isotropic Wavelet Transform on 
the sphere ( Starck et al. 2006a), and, for each wavelet scale, we 
have computed its second-, third- and fourth-order moment. We 
used 1 1 scales with central multipole values of £o = 9000, 4500, 
2250, 1125, 562, 282, 141, 71, 35, 18. For each of these maps, 
we computed the variance cr^ = (k^), the normalized skewness 
S = {K^)/cr^, and the normalized kurtosis K = {K^^y/cr^. Results 
are plotted in Fig. [3] as solid lines of various colors. Error bars 
were estimated approximately by computing each moment on 
the 12 Healpix base pixels independently and evaluating the vari- 
ance in the 12 results. A more appropriate strategy would have 
been to perform several, independent, 70 billion particle runs, 
which is currently impossible for us to do. We can see that the 
variance in the signal steadily increases for higher and higher 
multipoles, and saturates at a fraction of 10""^, corresponding to 
the value predicted from nonlinear gravitational clustering for 
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Fig. 5. Reconstruction of convergence maps with our 2 filtering techniques. The top panels show the 2.5^ x 2.5^ square map corre- 
sponding to first zoom sequence of Fig. [T] The bottom panels are subset of the corresponding top images with linear size 45'. From 
left to right, we show the original signal, the noisy image, the Wiener- filtered image and the the MRLens-filtered image. 



£ > 6000. The variance for each wavelet plane can be consid- 
ered to be a band power estimate of the angular power spectrum, 
as can be verified using Fig.lH In the same figure, we have also 
plotted for comparison the linear power spectrum, to highlight 
the scale below which nonlinear clustering contributes signifi- 
cantly, i^e^^_for_^_>_750_orequivalently < 15, as first pointed 
out by I Jain & SeljakI (1 19971) . Skewness and kurtosis are more di- 
rect estimators of the signal non-Gaussianity. Departures from 
Gaussianity occur around £ ^ 750 - 1500, where both statis- 
tics cross unity. Due to the large dynamical range of the Horizon 
simulation, we computed a map spanning two decades in angular 
scales in the linear, Gaussian regime and two additional decades 
in angular scales into the nonlinear, non-Gaussian regime. 

It is clear from Fig. [3] that at small £, the skewness and the 
kurtosis of the map are strongly aff'ected by cosmic variance. The 
statistics of the convergence field cannot be measured in practice 
over the whole sky because of sky cuts imposed by the pres- 
ence of saturated stars and by absorption in the Galactic plane. 
We estimated the impact of this sky cut on the accuracy of our 
multi-resolution statistical analysis. We computed the expected 
number of bright stars that would saturate CCD cameras typ- 
ically employed in wide-field survey (B-magnitude< 20). We 
then removed from our analysis each pixel contaminated by at 
least 3 bright stars, based on a random Poissonian realization of 
brig ht stars in our Galaxy (according to the model presented in 
iBah call & Soneira (1980), Appendix B). We obtain a mask with 
40% of the sky removed, corresponding roughly to a ±20^ cut 
around the Galactic plane (see Fig. O. The resulting statistics 
are overplotted as dotted lines in Fig.O The transition scale, for 
which the departure from Gaussianity is significant, can still be 
estimated reliabily around £ ^ 750 - 1500. We concluded that 
the cosmic variance of the cut sky aff'ects high-order moments 
only below £ ^ 200. 



5. Map-making using multi-resolution filtering 

The full-sky simulated convergence maps described above can 
be used to analyze and compare de-noising (or map-making) 
methods on the sphere. For this purpose, we considered a purely 
white instrumental noise, typical of the next generation all-sky 
surveys, and a root mean square per pixel of area Ap given 
by (Tn = 0.3/ ^rigaiAp for rigai = 40 background galaxies per 
arcmin^. Recovering the most accurate convergence map from 
noisy data will be an important step in future surveys. This re- 
constructed map can be used to construct a mass selected halo 
catalog, measure its statistical properties and constrain cosmo- 
logical parameters, and be compared directly with other cluster 
catalogues compiled with other techniques (X-ray, galaxy counts 
or SZ). We restrict ourselves to the full-sky denoising of a con- 
vergence map already reconstructed from the shear derived from 
galaxy ellipticities. In the present work, we do not address filter- 
ing in the presence of a cut-sky, such as the one shown in Fig.[2l 
Promising methods b ased on "impaintin g" have been developed 
in the CMB context ((Abrial et al.ll2008l) . and also weak-lensing 



applications (iPires et al.ll2008l) : these replace missing data with 
an artificial signal and allow us to optimize the results we ob- 
tained with filtering methods for a full-sky analysis. 

A straightforward filtering method is the Wiener filtering 
scheme, which is optimal for Gaussian random fields, and is ex- 
pected to operate here eff'ectively on large scales. Defining S{ sls 
the power spectrum of the input signal (see Fig. ^ and A^^ the 
power spectrum of the noise, this method involves convolving 
the noisy map by the Wiener filter defined as W{ = S{/(S{-\-N{). 
The results of the Wiener filtering approach are shown in Fig. [5] 
Comparing with the input signal map, we conclude that, al- 
though the agreement is satisfactory on large scales, the dense 
clumps clearly visible in the image are poorly recovered because 
they have been convolved too significantly. 
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Fig. 4. Angular power spectrum of the simulated convergence 
map (black solid line), compared to a fit based on the Smith et 
al. (2001) analytical model with error bars corresponding to our 
noise model (pink area). Also shown is the prediction from linear 
theory (pink dashed curve). The noise power spectrum is plotted 
as the dashed black line. The green solid line is the power spec- 
trum of the error map obtained with the Wiener filter method, 
while the blue solid line are that for the MRLens method. 



A dedicated weak-lensing wave let-restoration meth od, 
called MRLens, has been developed (IStarck et all l2006bb . It 
can be considered to be an extension of the Maximum Entropy 
Method (MEM) that provides diff'erent types of information. In 
MRLens, the entropy constraint is not applied to the pixels of 
the solution, but rather its wavelet coefficients. This allows us to 
take into account more efl&ciently the multi- scale behavior of the 
information. MRLens was, however, designed for weak-lensing 
maps of smaller surface area on the sky, for which the non- 
Gaussian signal is stronger. MRLens was extended here to the 
sphere by considering independently each of the 12 Healpix base 
pixels covering the sphere as 12 independent Car tesian maps, 
on whi ch we applied the MRLens algorithm of IStarcket aP 
(l2006 b!). Full- sky denoising performed with MRLens is shown 
in Fig. [5] It performs more efficiently than the Wiener methods 
on small scales, with dense clumps more accurately estimated, 
but less efficient than the Wiener method on large scale when 
recovering low frequency waves in the map. We also computed 
the angular power spectrum of the error map (see Fig.© in both 
cases (Wiener and MRLens). We can see that Wiener filtering 
outperforms MRLens on large scales. Interestingly, the MRLens 
errors decrease significantly above the transition scale we iden- 
tified in the last section around i ^ 1000 (see Fig. [4]). 

To compare both methods more quantitatively, we computed 
the skewness and kurtosis of both reconstructed maps. Results 
are shown in Fig.[6l We note that using map-making algorithms 
to recover the skewness and kurtosis of the true signal is not at all 
the optimal strategy: maximum likelihood estimators are more 
appropriate. We used high-order statistics here only to compare 
the relative merits of each method. It is striking to observe in 



Fig. 6. Skewness (blue lines) and kurtosis (red lines) for the orig- 
inal convergence map (solid lines with error bars), compared to 
the same high-order statistics for the Wiener reconstructed map 
(dotted lines) and the MRLens reconstructed map (dashed lines). 
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Fig. 7. Histogram of the residual maps for Wiener and MRLens 
filtering. 



Fig.[6]that the Wiener reconstructed map strongly underestimate 
the skewness and the kurtosis at small scale. This confirms quan- 
titatively what was already visible in the maps (Fig. [5]), namely 
that the Wiener method strongly suppresses high peaks in the 
map, aff'ecting the tail of the probability distribution function. 
On the other hand, the MRLens reconstructed map has a sig- 
nificantly higher skewness and kurtosis than the original map: 
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this wavelet-based method is only efficient in recovering high 
peaks in the signal, affecting the reconstructed probability den- 
sity function in the opposite direction. 

We now use the Probability Density Function (PDF) of the 
residual maps to compare each method (see Fig. [7]). We confirm 
our visual impression from Fig. [5] that MRLens performs more 
efficiently than the Wiener method in recovering the high conver- 
gence, nonlinear features in the map. The positive high residual 
tail is reduced significantly by MRLens, as well as the dozen 
of strongly outlying pixels in the Wiener filterer map around 
K ^ 0.35 (see Fig. [7]). MRLens, however, performs poorly for 
small values of the convergence (a: ^ ±0.05), for which the 
PDF is well approximated by a Gaussian, an optimal situation 
for Wiener filtering. 

The present analysis, based on using both the power spec- 
trum of the residual maps and the high-order moments of the 
recovered map, strongly suggests that new methods should be 
developed using an hybrid, multi-resolution formulation; for in- 
stance, using spherical harmonics on large scales, while utiliz- 
ing wavelets coefficients on small scales. The methodology of 
this combined approac h could be ba s ed on t he idea of Combined 
Filtering introduced by lStarck et al.l (l2006ah . 



6. Conclusion 

Using the 70 billion particles of the Horizon N-body Simulation, 
we have computed for the first time a realistic full- sky conver- 
gence map with a pixel resolution of A^ ^ 0.74 arcmin^. We 
have analyzed the resulting map using multi-resolution statistics 
(variance, skewness, and kurtosis) and angular power- spectrum 
analysis. We have shown that this simulated map spans 4 decades 
of useful signal in angular scale, with 2 decades within the lin- 
ear, Gaussian regime and 2 decades well into the nonlinear, non- 
Gaussian regime. We have shown that, when considering a real- 
istic sky cut, we can reliabily estimate high-order moments of 
the map above I ^ 200. Using even higher resolution maps, 
angular scales smaller than 6 ^ \' could be explored in future 
works, although the mass d itribution on the se scales might be 
affected by baryons physics (iJing et al.ll2006h . so that the present 
map might already cover all cosmologically relevant scales. 

As a first step towards a realistic map-making procedure, 
we have tested two de-noising schemes on a simplified fictitious 
dataset derived from the f ull- sky map, namely Wiener filtering 
and the MRLens method (IStarck et al.ll2006bl) . We have shown 
quantitatively that Wiener filtering is the most effective method 
on large scales, although some signal is lost on small scales. 
MRLens performs more effectively on small scales and recovers 
the dense clumps associated with dark matter halos, but deals 
less accurately with low frequency waves in the map. Hence, 
this work demonstrates the need for hybrid multi-resolution ap- 
proach, e.g., by combining spherical harmonics and wavelet co- 
efficients. The present analysis will be extended in future work to 
map-making algorithms dealing directly with galaxy shears. The 
simulated convergence map may prove to be an effective tool for 
the design of new map-making methods and for the preparation 
of the next generation weak-lensing survey^. 
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Appendix A: Computing the convergence maps 
from simulations 

We first recall how to compute the convergence in the Born ap- 
proximation, and then present our new ray-tracing scheme. 



A.1. Born approximation 

We start by the formula relating the convergence to the density 
contrast: 



3 r 



dz D{z)D{z.Zs) ^ ..c 

-S(-—D(z)n,z), 
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which is valid for sources at a single redshift Zs, and D(z) = 
^X(z) is the dimensionless, comoving, radial coordinate (dD = 
dz/E(z)). We now rewrite this formula in a form that is more 
suited to integration over redshift slices in a simulation: 
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is a slice-related weight, and the integral over the density con- 
trast reads 
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where 

Spi.{Zb) = ^—-^D\zb) 
lypix /7q 

is the comoving surface of the spherical pixel. Interpreting all 
together, we obtain the following formula for the convergence 
map (omitting the A;^/, term that corresponds to a constant term): 
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This is the equation used to derive the convergence map in the 
Born approximation. 

A.2. Ray-traclg using multiple planes 

We discuss here the formulae needed for the multi-plane com- 
putations, where we consider the lensing by a number of thin 
lenses located at {zb}- We define 
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^{Zb,0) = KfacOj(Zb)- 



Npart(0,Zb) 
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with 



M{Zb) 



( r dz n(z) \ 



Jaz. E(z)) 



To follow the light rays, we are interested in computing the an- 
gular displacement field for each ray / due to a slice at zt- We 
then define 



a' 



-2VA- 



(azb)))m. 



(A3) 



where the gradient and Laplacien are computed using angular 
covariant derivatives on the (unit) sphere, and Ot is the current 
direction of light ray / when it is incident on the slice b. Now, we 
start from light rays that are back-propagated from the observer 
at z=0 towards the source (here at z=l). We denote by {0^} the 
location of the Healpix centres, which corresponds to the initial 
directions of the back-propagated rays emanating from the ob- 
server. The tangent vectors to each light ray will be modified by 
the deflection field at each lens plane, defined by Eq. IA.3I Then, 
computing the displacement of the rays at slice b reads 

a^ = {-2VA-\azb))(0^)- 

We then update the direction )0^ of the rays according to the fol- 
lowing rotation, "R: 



;8f = ^(«fxaf,||af||);8f- 



(A.4) 



where J3^ = n^^ (light rays emanate from the observer, thus in 
a direction perpendicular to the first slice), and fi^is the vec- 
tor normal to slice b at the intersection of light-ray / on slice 
b. Equation (IA.4I) can be simplified by noting that a is expressed 
naturally in the local (eg ^0) basis of the tangent plane at position 

(o,cpy. 
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a = ||a||(cos^^6' + sin^^0) . 

After calculating the new value of )0, one needs to compute the 
intersection of the light rays with the next shell. We call x^ the 
Cartesian position of the intersection of light ray / with slice b, 
then the next intersection will be given by the solution for A of 
the system: 



A^^2A(x^-fi^)^Rl 



>+i 



-^b+i 



= X) 







A>0, 



assuming that fi remains strictly unitary, and Rb is the comov- 
ing radius of slice b. Once jc^^^ is known, it is easy to compute 
the new 6^'^^ positions. The contributions to k are then calculated 
following Eq. lA.ll but where the slice contributions are interpo- 
lated at the displaced positions: 



m ^ ^n ^p- / ^0 f V(sim) ^ 

2 4n \ c / Npartisim) ^ 



Npart(e^,Zb) 

mzb) 



We note that ^ may fall into different pixels as a function of the 
slice b. 



